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Abstract. We study the evolution of a system of free fermions in one dimension 
under the simultaneous effects of coherent tunneling and stochastic Markovian noise. 
We identify a class of noise terms where a hierarchy of decoupled equations for the 
correlation functions emerges. In the special case of incoherent, nearest-neighbour 
hopping the equation for the two-point functions is solved explicitly. The Green's 
function for the particle density is obtained analytically and a timescale is identified 
where a crossover from ballistic to diffusive behaviour takes place. The result can 
be interpreted as a competition between the two types of conduction channels where 
diffusion dominates on large timescales. 



1. Introduction 

Transport properties of one-dimensional quantum systems sliow intriguing features and 
have been the topic of intensive research. Despite the efforts, some important aspects of 
the problem still lack a conclusive answer. The main open question is whether ballistic 
transport, characteristic of a number of integrable quantum systems at zero temperature, 
could survive thermal noise or rather a transition to diffusion is inevitable. Paradigmatic 
integrable model systems include weakly interacting fermions or, in the context of spin 
systems, the XXZ model where a number of analytical methods are available [Tj. In 
spite of the long-standing conjecture that integrability protects the ballistic features of 
transport at any finite temperature [2] , recent calculations have pointed out the existence 
of a dominant diffusive transport channel in the gapless phase of the XXZ model at half 
filling [3l H] . The above results, along with a large number of different numerical and 
analytical works [5], are based on linear response theory and the calculation of the 
conductivity through Kubo's formula. 

A completely different approach to the transport problem is possible in the 
framework of open quantum systems [6]. Under certain assumptions on the coupling 
to the environment, the time evolution can be cast in the form of a quantum master 
equation for the system density matrix. In the context of spin chains this approach 
was initiated in [7] by modeling the coupling of the chain at both ends to heat baths of 
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different temperatures. In the Markovian approximation, an analytical treatment of the 
problem was presented [8] by diagonalizing the time evolution operator for spin chains 
that can be mapped into a quadratic fermionic system. In the steady state of the XX 
chain (respectively free fermions) a fiat magnetization profile emerges [9], suggesting a 
ballistic transport. This has been recently supported by a calculation of a lower bound 
on the Drude weight for the critical XXZ chain [TO] . 

The above approach is restricted to incoherent processes which either create or 
annihilate a particle. Furthermore, most of the cases considered so far were restricted to 
boundary driving. An interesting example, where particle-conserving dephasing noise 
is present in the bulk; of an XX chain was presented in [TT]. The steady state was 
calculated analytically in a perturbation series with respect to the boundary driving 
and some exact results show a diffusive, linear magnetization profile for any finite value 
of the bulk noise. Similar results were found numerically for the XXZ chain [12] as well 
as in the case when the noise term describes stochastic hopping [13] . 

It must be emphasized, that most of the above mentioned examples considered only 
steady-state properties and the few existing results on the time evolution are limited 
to numerical methods [HI ITS] . On the other hand, there exists a remarkable exact 
solution for the complete spectrum of the time evolution operator of a quantum diffusion 
problem [161 [T7]. This is, however, restricted to a single fermion moving on a chain with 
coherent tunneling and subject to dephasing noise. Therefore, it would be favorable 
to find examples where time evolution under particle-conserving stochastic noise can 
be exactly tackled in the context of a genuine many-particle problem. This is further 
motivated by the fact, that such models have been recently suggested for the description 
of energy transfer in photosynthetic complexes and biomolecules [HI [19]. 

In the present paper we discuss the evolution equations of free fermions when the 
coherent tunneling is supplemented by thermally activated, stochastic hopping. This 
choice is favored since in the infinite temperature limit it reproduces the well-known 
symmetric simple exclusion process [20]. An investigation of the evolution operator 
shows a remarkable hierarchy property of the correlation functions which enables us to 
derive the exact equations of motion for the two-point correlations. These equations 
turn out to be very similar to the case of the single particle quantum diffusion [16] 
and can be explicitly solved. A further analysis of the resulting expressions for the 
time-dependent density profile yields an exact analytical form of the Green's function. 
This formula is then used to identify a well-defined timescale which separates a ballistic 
transport regime for short times from the diffusive behaviour in the long time limit. 
The effective diffusion constant diverges in the limit of vanishing noise, signaling the 
transition to the pure ballistic regime. 

In Section [2] we investigate the master equation for the density matrix and introduce 
a class of stochastic processes where the structure considerably simplifies and a hierarchy 
of decoupled equations emerges. The simple case where the stochastic terms describe 
symmetric nearest-neighbour exclusion is introduced in Sec. Ejfollowed by the derivation 
of the master equations for the two-point correlations. Section H] is devoted to the 
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analytical solution of these equations. In Section E] we give an explicit analytical form 
of the Green's function for the density profile and compare it with numerical results. 
Our findings are discussed in the last Section [6] and some details of the calculations are 
presented in two Appendices. 

2. The master equation 

The coherent evolution of closed quantum systems is described by a unitary time 
evolution operator. However, in any realistic situation one has an inevitable coupling 
to the environment which introduces incoherent effects. Under certain assumptions, the 
description of these open quantum systems is possible in terms of a quantum master 
equation involving the density matrix of the system [6]. In many examples the dynamics 
of the environment has a much shorter time scale and it is reasonable to make the 
Markovian approximation. Then the time evolution of the system density matrix p is, 
in general, given by a master equation of the Lindblad form [21] 



where the symbol C refers to the Liouvillian. The first term in the equation corresponds 
to a coherent time evolution according to Hamiltonian H, while the operators La 
describe different sources of stochastic noise. 

We will focus on fermionic quantum systems defined on a one-dimensional chain 
of length N. The degrees of freedom are described by the creation and annihilation 
operators and aj„ possessing canonical anticommutational relations {0^,0]^} = 5m,n 
and {am, an} = for m, n = 1, . . . , A^. We will be interested in model systems where 
the coherent evolution is given by a free fermionic Hamiltonian that is quadratic in 
the creation and annihilation operators. However, without further assumptions on the 
Lindblad operators La, solving the master equation ([T]) turns out to be a very difficult 
problem. 

Recently a simple integrable example was given by Prosen [8]. Here, the Lindblad 
operators were taken to be an arbitrary linear combination of the fermionic operators 
and in turn describe processes involving the creation or the loss of a particle. The 
integrability of the master equation relies on the fact that the operators La appear 
quadratically in ([1]) and induce a Gaussian time evolution operator. This is, however, 
no longer true if one considers stochastic processes (e.g. simple exclusion) that conserve 
the number of particles. Then La involve quadratic terms in the fermi operators and 
the Liouvillian is not any more diagonalizable by a canonical transformation. However, 
as it will be shown below, some particular choice for the incoherent terms leads to a 
considerable simplification of the problem. 

We will follow the formulation of Ref. [H] and start by introducing 2N Majorana 
fermions with the definition 




a 
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'm 1 
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that are Hermitian and satisfying the relations {c^, q} = 2Sk^i for all fc, / = 1, . . . , 2N. In 
terms of these operators the most general quadratic Hamiltonian and Lindblad operators 
read 

^ 2N ^ 2N 

-f^ = ^ ^ HklCkQ , La = - ^ La^klCkCl (3) 

k,l=l k,l=l 

where Hik = —Hki is required by hermiticity of the Hamiltonian. Note, that this choice 
of H and La corresponds to coherent and incoherent processes that either conserve the 
particle number or involve pair creation and annihilation. 

It turns out to be more convenient to work with observables instead of density 
operators. For this purpose one can introduce a convenient basis and define the ordered 
strings of Majorana operators 

r, = cr...c^?7, z/, €{0,1} (4) 

where z/ = (z/i, . . . , z/2Ar) denotes the vector of the occupation numbers z/j indicating 
whether the corresponding Majorana operator q is present in the string Ty^. These 
objects encode the different n-point correlation functions where the order is given by 
n = z/j. It is also useful to define superoperators acting on these strings that create 
or annihilate a Majorana operator at position j as 

CjT^ = 5i^u,T^jT^,, c\Ti, = So^u.T^jTi,/, ul = \ ^ ^. .{ (5) 

where the sign factor tt^ = exp {it: X]i=i ensures that canonical anticommutational 

relations {Q,cj} = and {ci,Cj} = are satisfied. Note, that these superoperators 
change the order of the correlation function from n to n ± 1. The overall number of 
independent, ordered correlation functions is 2^^ which exactly equals the number of 
components in the density matrix p. This indicates, that solving the master equation 
([1]) for p is equivalent with solving the complete set of master equations for the operators 
Vtj_ which are generated by the adjoint Liouvillian . 

We are now ready to state the main result of this section. If the Lindblad operators 
satisfy L]^ = La for all a the Liouvillian takes the following simple form 

2N ^ 2N 

= -Yl Hkiclci + 2 5Z (^) 

k,l=l i,j,k,l=l a 

where Hki denote the elements of the matrix H = H + | J2a L^La- It is important to 
stress that the operator ([6]) conserves the number of Majorana fermions and therefore 
the length of the strings Ty^. In other words, a hierarchy of the n-point correlation 
functions emerges since the time evolution under does not mix strings of different 
length. Note, that a similar hierarchy property was pointed out for the steady state of 
the driven XX chain with dephasing [22] , which is now seen to generalize to the complete 
time evolution. For the derivation of ([6]) we refer to Appendix A. 
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Before fixing our model, it is worth investigating the general structure of Eq. 
The first term contains the Hamiltonian evolution modified by a damping term which is 
in turn responsible for the decay of the n-point functions. This quadratic term generates 
the Gaussian part of the time evolution and has exactly the same form as in case of 
linear Lindblad operators [8] . However, one has an additional "interaction term" which 
is of fourth order in the superoperators and in general implies that time evolution is 
non-Gaussian. Although the Liouvillian cannot be diagonalized in general as in [8], 
some simple choice of the stochastic terms could further simplify ([6]) and eventually 
lead to a tractable problem. 

3. The symmetric quantum exclusion process 

After setting up the general formalism, we proceed to specifying the concrete model. Our 
focus is to define a dynamics which interpolates between the extreme cases of quantum 
coherent tunneling, described by the tight-binding Hamiltonian, and classical stochastic 
hopping, described by the symmetric simple exclusion process [20]. Such a model has 
been introduced recently [13] and investigated numerically from the perspective of steady 
state properties. In order to get a nontrivial steady state, boundary injection and 
ejection was introduced into the dynamics. Since here we are interested in the complete 
time- dependent solution, it turns out to be more convenient to take periodic boundary 
conditions. The geometry and the update rules are sketched in Figure [H 




Figure 1. Symmetric quantum exclusion process on a ring of size N. Particles can 
move to unoccupied nearest neighbour sites by coherent tunneling (curled arrows) or 
stochastic hopping (simple arrows) with corresponding rates A and 7. 

The coherent evolution is described by the Hamiltonian 

N 

H = -X ^{a]aj+i + a]^^aj) (7) 
i=i 

where A is the tunneling rate. The incoherent hopping is generated by the Lindblad 
operators 

^Lj = , Lrj = ^a\^^aj (8) 

where Lj^j and L^j with j = 1, . . . , describe hopping to the left and right neighbour 
sites with the same rate 7. Although these Lindblad operators are not Hermitian, one 
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can introduce new ones with a simple unitary transformation 

1 i 
L2j^i = {Llj + Ljij) , L2j = {Llj - Lrj) (9) 

which now satisfy the conditions L„ = Lj^ for every a = 1, . . . , 2N. Since the master 
equation ([1]) is invariant under such unitary transformations, one can now apply the 
results of the previous section. Note, that this step requires symmetric hopping. 

The matrix elements of H in the Majorana basis can be written in a block matrix 
notation as 

Hmn ^ ^ iSm,n-l + Sm,n+l) (10) 

where the indices take the values m,n = 1, . . . ,N. Similarly, one has 

;)^(^.A...-^™...M (11) 

The last matrix is diagonal and describes a pure damping term. Moreover, all 
the other matrices have a simple block-tridiagonal form, however some non-diagonal 
entries appear. Therefore, we perform the following canonical transformation of the 
superoperators 

flm = (C2m-1 ~ 'i'C2m) , = (C2m-1 + ^C2m) (12) 

which will diagonalize the 2x2 matrices in H and i^2j-i with eigenvalues ±i. Note, 
that the effective action of the superoperators dm and bm is to remove the fermionic 
operators and a]^ from a string. 

Substituting into (|6]) a simple calculation yields the Liouvillian of the symmetric 
quantum exclusion process in the form = C^a + + i^'^ ab where 





m 




m 




m 



(13) 



-iblnbm+l + bln+ib 



it 



The operator a is therefore formally equivalent to the one describing the motion of 
a ID fermionic system in a chemical potential and with repulsive nearest neighbour 
interactions. Note, however, that one has an imaginary hopping amplitude i\. The 
term C\ describes a second species of fermions where the only difference is in the 
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hopping amplitude —iX. The last term translates into a special interaction between the 
two types of fermions, giving a penalty term for simultaneous hopping. 

The time evolution of the correlators of the quantum exclusion process is 
obtained by systematically applying (fT3|) to the various strings of operators, which 
is demonstrated below on the simplest examples. 



3.1. One-point correlations 

The simplest string consists of only one fermion operator. The contributions of the 
interaction terms then vanish and one has 

= C^aicim) = (a^^i + am+i) - 'yam (14) 

while the equation for follows simply by hermitian conjugation. Using the Fourier 
transformed operators a^, this equation can be solved as 

ttg = e^*ag(0), e = 2iXcosq — -yt (15) 

where q = 27m/N for n = 1,...,N. Since most of the physically interesting initial 
states will give a vanishing expectation value (ag(0)), we move forward to analyze the 
first nontrivial correlation functions. 



3.2. Two-point correlations 

The most basic physical quantities, such as the particle density, are encoded in the 
two-point correlations. The independent correlators are chosen as aj^an for m <n and 
al^a\ for m < n, while all the other combinations are related by complex conjugation 
and the commutational relations. Instead of dealing with operators, one can already 
take the expectation values with respect to some arbitrary initial state. The equation 
for Gm,n = {(^In'^n) then reads 

1 , GjYi.n ^A {^Gm—l.n GjYi+l.n G^.n—l G^.n+l) '^iG jyi.n „n 

at (16) 

where the contribution in the second line is generated by ab- Similarly, the evolution 
of the pair-creation expectation values Fm,n = {'An'Ai) given by 

-Fm,n = -lX iFm-l,n + + Fm,n-1 + Fm,n+l) ' 2^Fm,n ^^^^ 

~l~ 'y^m,n—lFra,m+l ^^m,n+lFin— 1,171 

Since both ( !T6|l and ( !T7|1 are linear, first order differential equations, they can be 
formulated as a matrix eigenvalue problem. The interaction term of the Liouvillian has 
only a localized contribution proportional to 6m,n and Sm,n±i, respectively, and therefore 
the equations translate into a simple potential scattering problem. 
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4. Solution of the scattering problem 

The equation (fT6l) for the particle-hole correlators has a similar form to the one appearing 
in [L61 as the equation of motion for a single electron with dephasing noise. The solution 
of the scattering problem is presented in detail for Gm,n along the lines of Ref. [16] and 
we also briefly comment on the analogous treatment for Fm^n- 



4.1. Particle-hole correlations 

The linearity of ( IT6l) allows one to separate the time dependent part as Gm,n(t) = Gm,ne!^^ 
with the energy eigenvalue e. Translational invariance implies, that the eigenvectors 
obey Gm,n = e*''"'G'o,n-m where the allowed values of the wavenumber are qj = 
with j = 1, . . . ,N. Introducing the relative coordinate I = n — m and substituting 
Go.z = i^^e^^^gi we get the simpler equation 

q 

egi = 2i\sm-{gi_i + gi+i) - 2-fgi + 2-f cos q giSi^ (18) 

where < / < — 1 and the boundary conditions imply g^ = i^ e~^^^ go. One 
recognizes, that the structure of Eq. (ITS!) for the amplitude gi is the same as the one 
describing the motion of a free particle with a potential scattering center localized at 
/ = 0. This is solved by using the ansatz 

gi = Ae'^^ + Be-'^^ (19) 

Substituting into (fT8|) one obtains A^— 2 identical equations that fix the energy eigenvalue 
as 

Ege = AiX sin | cos ^ - 27 (20) 

where q and 6 are explicitly used to index the eigenvalues. The remaining two equations 
for / = and / = — 1 determine the quantization of the second quantum number 6 
and fix the value B/A of the scattering phase. The solution of these equations follows 
exactly along the lines of Ref. [H] and is summarized in Appendix B. In turn, one finds 
that the eigenmodes can be written in the following Bethe ansatz form 

G^J^„ = zTzl^ + S,^{-Z2r{-z,r (21) 
where S21 = B/A and we used the notation 

= g%/2-e+V2)^ _ g%/2+e-V2) _ (22) 

The allowed wavenumbers 6 for finite A^ can only be determined numerically by 
solving Eq. (IB.ip . The resulting zi and Z2 are shown in Fig. |2]for A^ = 70. In spite of 



the nonvanishing imaginary parts of 6, one clearly sees a condensation around the unit 
circle, corresponding to the eigenvalue families 61 and 611. However, the third family 
6d leads to a separate branch of solutions Zi and Z2 on the complex plane, lying either 
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Figure 2. Allowed values of the parameters z\ (red dots) and (blue dots) on the 
complex plane for N = 70. The tunneling rate is A = 1 and the hopping rate takes 
the values 7 = 0.1 (left), 7 = 1 (middle) and 7 = 8 (right). The zi solutions of the 
diffusive branch lie further outside of the unit circle and are not shown. 



inside or outside (not shown in figure) of tlie unit circle. Tliese solutions correspond to 
purely real energy eigenvalues 



£rf = 27jcos2g-4^sin2|- 27. (23) 



It is obvious from Fig. [2] that their number increases for larger values of the stochastic 
noise 7 and the solution 2:2 = in the center of the circle, corresponding to the steady 
state = 0, is always present. These eigenvalues are, in turn, responsible for the 
emergence of diffusion in the time evolution. It is intuitively seen by expanding (1231) 
around the steady state g = which gives 

€d ~ -Dq^, D = ^ + — (24) 

7 

and therefore yields the usual form of a diffusive dispersion relation with diffusion 
coefficient D. For the detailed discussion of the evolution of the particle density we 
refer to the next section. 

In order to obtain the general solution matching to a given initial condition one 
has to first construct the left eigenvectors and prove the orthonormality relation. One 
starts by noticing that the spectrum in Eq. ( I20l) is degenerate for the pair of solutions 
g — )■ 277 — g while the complex phases transform as Zi — )■ (—-22)"^ and Z2 — )■ {—zi)~^ 
under this change. Therefore, we write the left eigenvector as 

G'm'n = ^r'"^2-" + S,:{-Z,r-^{-Z^)-- (25) 

In order to obtain an orthonormal system one has to extend the eigenvectors to 
index pairs where m > n. Since the solutions of equation (1181) are either even 
g^i = gi or odd g-i = —gi, the right eigenvectors must satisfy the symmetry property 
G'^f^ = ±(— 1)™+"G'^^„ and the same holds for the left eigenvectors. Note, that the 
antisymmetric eigenvectors, insensitive to the value of the scattering potential, are 
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exactly the ones in the family On (see Appendix B). Introducing the notation (q', 6'\ for 
the left and \q,0) for the right eigenvectors, a lengthy but simple calculation shows 

{q', e'\q, e) = J2 GtiG^l^ = 5,,,.5e,e'N,e (26) 

m,n 

where the sum goes over all the indices and the complex normalization factor reads 

iV,, = 2iV(iV + .A), A = -7^^^^ (27) 

1 — sm d 

where the parameter /3 is defined in (]B.2|) . Finally, the complete solution can be written 
in the form 

m)) = Y.^,e\q,9)e^^^\ c,e = N;,'{q,e\G{0)) (28) 
q,e 

where the constants Cqo are determined through the initial condition |G(0)). 

To conclude this section, we show that the inverse of the complex normalization 
constant A^^^ can be interpreted as the density of states in the wavenumber space. First 
we note, that the solutions 6 have a nonvanishing imaginary part which depends on the 
real part and thus also changes when going from one solution to the next. Using the 
notation of Appendix B one has d6 = ddo + id6 and by differentiating ( IB. 51) one finds 
d^^ = 2-n[N + iA)^^. Therefore one has A"^^ = dgd0/87r^ which, up to a factor, indeed 
corresponds to the level density. 



4-.2. Particle-particle correlations 

For completeness, we also remark on the solution of the pair-creation probabilities 
Fm,n(t). In principle, one has to follow the same steps as in case of the particle- 
hole correlations. The time dependence is again of the form Fm,n(t) = Fm^n^-^^ and 
the complex parameters now have to be defined as zi = e*'-'^/^"^-' and Z2 = e^^'^^'^~^^\ 
respectively. It is easy to check that the ansatz 

Fl'n = ^T^2 + S21ZTZ^ (29) 

solves (I2T]) if the energy eigenvalue and the scattering phase satisfy 

~ q ~ Z\Z2 + 1 — i\z2 

Sqe = -4:iX cos -cos 6 -2^, 821 = --^-: ^tt- . (30) 

2 ^12:2 + 1 - ij^zi 

The allowed values of the wavenumber 6 is again determined by a corresponding Bethe 
equation which we do not discuss in detail. Finally, one should remark that the 
scattering phase (!30ll has a very similar structure to the one appearing in the Bethe 
ansatz solution of the simple exclusion process [23] . 



5. Green function of the density 

In this section we will focus on the time evolution of the particle density Gm,m{t)- The 
essential quantity to be determined is the Green's function Q{m — k,t), describing the 
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evolution of a single particle localized at site k for t = 0. The initial condition then 
reads 

Gm,n{0) = Sm,kSn,k (31) 

and the scalar products with the left eigenvectors yield 

(g,^^|G(0)) = (l + 52ll)e-'?^ (32) 

Therefore, the family On where 5*21 = —1 does not enter the solution. Using the form 
flB.4p of the scattering phase the evolution of the density reads 



a(m - *, = ^ ^"t^S^ (33) 

q,e 

where the sum has to be taken over the families Oj and Od- As remarked at the end 
of section 14.11 the factor A*"^^ naturally translates into the density of states in the 
N ^ oo limit. The diffusive eigenvalues satisfy P sin 6 a 1 but at the same time have 
a vanishing density N~q^ — )■ in the 6'-space thus the product of the first two factors in 
fl33|) is finite and can be evaluated using fl27|) . Hence, the Green's function is a sum of 
two separate contributions 

Jqx 1 ^egt 



'|/3|<i 27r ^1 - /32 

r d0 (S^sin^e 



lo 271 Jo n 1- 13^ sin^ 9 
where we introduced x = m — k. The first integral goes only over q values satisfying 
< 1 and the ^-integral has a pole for \/3\ > 1, corresponding to the branch cut in the 
9i eigenvalues. This suggests writing 

c(M)^r?e(,,oe"-. efe,*) 4 (35) 



27,-"' ' ' [ g'{q,t), if |/3|>1 

where the Fourier transformed Green's function has to be defined piecewise. It must 
be emphasized, however, that Q{q,t) turns out to be a smooth and continuous function 
and its different forms on both sides of = 1 are connected by analytical continuation. 

In order to obtain (g, t) , it is first useful to introduce the parameter a = cos q, 
in terms of which one has Ed = ol\/\ — (5'^ — 2'yt and Ego = aP cos6 — 2'yt. For < 1 the 
denominator of the 6'-integrand can be expanded in a geometric series. The resulting 
integrals can be carried out |2lj and yield 

Jn{oiP){2n 



e~'^' . (36) 



The above functional form is obviously only valid for /3 < 1. In spite of the apparent 
singularity of the first term for /3 — )■ 1, the sum also becomes divergent for these values 
and, since it enters with the minus sign, regularizes the expression. In order to see this 
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one can expand the exponential in terms of its variables a and /3 and further use a series 
representation of the Bessel function [23] 

f^Vy:M);,f^V' (3.) 



k=0 



k\{n + k)\ 



which then yields Q^{q-, t) as a double infinite sum. One can see that many terms cancel 
out and, on one hand, one obtains the first term of (1361) with the exponential replaced 
by a sinh function. Additionally, one has a regularized series of the same expression 
with cosh, where one of the sums is cut at a finite order. With a proper reordering of 
the latter terms and using again flTTI) one can rewrite it as a single sum involving Bessel 
functions. Finally, the term with the hyperbolic sine can be analytically continued to 
> 1 and one arrives to the following expression 



sm a 



n=0 



{2n - 1)! 



-2yt 



(3J 



Note, that the first term is just the pole contribution of the integral ( l34l) . Since the 
result (138|1 is obtained by analytic continuation rather then evaluating (134|) for > 1, 
the continuity of the function Q (g, t) is guaranteed. 

The final Fourier transform (l35l) which yields the Green's function in coordinate 
space in general has to be evaluated numerically. However, on some well defined 
timescales one finds simple approximations to G{q,t) and the integral can be carried 
out explicitly. 
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Figure 3. Fourier transform Q{q,t) of tlie Green's function for slrort (left) and long 
times (right) with parameter values A = 1 and 7 ~ 0.1. Dashed lines correspond to 
the approximations (|39p and (j41[) . respectively. Points represent the exact solution for 
finite iV = 70 by evaluating the sum over in 



Short time behaviour 

On short timescales compared to the stochastic hopping rate, 7^ ^ 1, one finds that for 
all wavenumbers satisfying the condition a <^ /3 only the n = term in the sum of Eq. 
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contributes significantly and one has 



sin(a/3) ^ , ^, 



(39) 



/3 

The left hand side of Figure [3] shows that (!39|) indeed gives a very good approximation 
of Q {q, t) on short timescales. However, as time increases the deviation becomes large 
around g ^ where < 1 and one cannot use the representation Q^{q,t). 

The arguments in fl39|) can be rewritten as a/? = 4 At sin |. By further assuming 
7 ^ A, the first term can be neglected since it is multiplied by an additional factor of 

~ 7/A. The remaining term can now be explicitly integrated as [21] 

g{x,t)^[U2\t)fe-^'". (40) 

The resulting Q{x, t) is simply the Green's function of the coherent case multiplied by an 
exponential damping factor. Therefore, coherent effects are washed out on a timescale 

Long time behaviour 

An other simple approximation can be obtained for large times 7^ ^ 1. Then the sum 
in (!38|l can be neglected unless a /3 where one has to use the other representation 
Q'^{q,t). Since in (15^ the n = term is missing from the sum, one has 

(41) 

The right hand side of Fig. 12] shows that for increasing times the oscillating tail of 
G{q,t) is suppressed and the approximation works well for the nonvanishing part at 
small wavenumbers. Finally, using the expansion fl24p of Ed and noticing that for 'yt ^ 1 
the denominator in fHTl) can be set equal to 1, one arrives to the simple result 



Therefore, the long time behaviour of the evolution is always governed by diffusion. 
Nevertheless, one sees from the expression of the diffusion constant ( 124]) that the speed 
of diffusion can be significantly enhanced and actually diverges as 7 — 0, signaling the 
transition to pure ballistic transport. 



Finally, we compare the exact numerical results of the density evolution obtained 
for a chain of finite length = 70 with the numerically evaluated integrals for Q{x,t). 
As shown in Fig. H] the data show excellent agreement and verifies the analytical result 
for the Green's function. On the left hand side a relatively small stochastic hopping 
rate 7 = 0.1 allows the ballistic features to survive for short times, but for larger times 
the crossover to a Gaussian profile can be observed. The discrepancy of the finite size 
results for t = 20 is a simple consequence of the ring geometry, since the two ends of 
the expanding wavefront collide. The right hand side shows results for 7 = 1, where 
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already the smallest time shown is in the diffusive regime. The difference between the 
speeds of spreading is clearly visible. 




Figure 4. Green's function of the density for different times and parameter values 
7 = 0.1 (left) and 7 = 1 (right). The tunneling rate is set to A = 1. Points represent 
the exact solution for finite iV = 70 by evaluating the sums in (|33|) . 

6. Discussion 

We have presented an exact analytical solution for the Green's function of the density 
profile in the symmetric quantum exclusion process. The derivation is based on a more 
general framework for the treatment of the quantum master equation for stochastic 
processes with quadratic and hermitian Lindblad generators. In our specific model the 
Liouvillian becomes analogous to an operator which describes a coupled system of two 
species of fermions with nearest-neighbour interactions. This form allowed an exact 
treatment of the one- and two-point correlation functions. 

For the complete solution of the problem one should also look at the higher order 
correlations. The form of the Liouvillian (fT3l) suggests that the solution should be 
available by Bethe ansatz. This would generate the n-point functions in a form similar 
to fl2T]) and fl29|) but with additional terms where all the different permutations of the 
phase factors are present. The corresponding scattering phases should then factorize into 
pair-terms, giving a factor Sji when the permutation exchanges a particle at position i 
with a hole at position j and a factor a Sji when two particle-operators are exchanged. 
Although this scheme looks feasible, a rigorous proof such as the one recently given for 
the simple exclusion process [25] would be desirable. 

It is important to stress that the short and long time behaviour of the density 
results from the asymptotic behaviour of the Green's function for large and small 
wavelengths, respectively. Therefore, it can be interpreted as a competition between 
ballistic and diffusive channels for the spectral weight in the propagator. This is exactly 
the same mechanism which was outlined in [4j by the discussion of the thermal spin-spin 
correlation function for the XXZ chain. Note, that in our case the large wavenumbers 
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are always exponentially suppressed for long times and therefore the diffusive transport 
channel dominates. 

With the solution for the Green's function in hand, one could further look at 
evolutions from some simple initial conditions, such as a step function in the density 
profile. This situation was investigated recently for the symmetric simple exclusion 
process |26] as well as in the purely coherent case [27t [28] . Apart from the crossover 
in the evolution of densities, one could look at the fluctuations of the particle number 
where an interesting crossover from y/t [26] to logt scaling [28] should emerge. 

Although the problem was formulated in the language of fermions, it is 
straightforward to generalize it to certain spin models, such as the XY chain, which 
can be transformed into quadratic fermionic Hamiltonians. However, one also has to 
fulfill the requirement that the stochastic terms translate into quadratic and hermitian 
Lindblad operators. The simplest example is the so-called dephasing noise, which is 
described by a term and is therefore included in this class. 

Another extension of the present work could be considering different boundary 
conditions. It would be interesting to check whether the problem remains solvable if 
one has, instead of a ring, a linear chain with boundary injection and ejection at the 
ends. Since an exact solution for the steady state has been recently presented in case of 
the XX model with dephasing noise [HI [22] , one could speculate whether our framework 
could generalize results considering the entire dynamics. 

The treatment of the two-point correlation functions in the present many-body 
problem greatly parallels to the single-particle quantum diffusion problem of Esposito 
and Gaspard [16]. The only difference is in the actual form of the function /3 which, 
however, does not change the qualitative picture and, in the thermodynamic limit, 
leads to diffusion on large timescales. We expect that this behaviour would not change 
by considering a stochastic process with longer range hopping, although the effective 
diffusion constant might increase. Therefore it remains a puzzling question, whether 
one could construct some more complicated stochastic processes which would eventually 
protect the ballistic features and lead to the onset of diffusion at a finite value of the 
corresponding rate. 

Finally, one should point out that the hermiticity of the Lindblad operators was 
found to be a sufficient condition for the existence of the correlation function hierarchy. 
It would be interesting to address the question whether it is also a necessary one. 
Furthermore, one could investigate whether the inclusion of some interaction term in 
the Hamiltonian would still leave the hierarchy unchanged. 
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Appendix A. Time evolution of the n-point correlations 

In this appendix we derive the formulas that are used to obtain the time evolution 
of the ra-point functions Fj, introduced in Eq. (jl]). Using the commutation relation 

{cfc, q} = 26k,i it is easy to show that 

[ckCi, r,] = 2 (clci - cjcfc) (A.l) 

where the superoperators and Ck that insert or remove a Majorana operator from the 
string Fj, were defined in ([5]). Together with the definition of the Hamiltonian in Eq. 
([3]) and Hm = —Hik this yields 

2N 

4oh(r^) = nH,r._] = -Y. HkidiciT^ (A.2) 

k,l=l 

and one obtains the coherent part of the evolution. Setting C = Ccoh + '^a, one has 
to calculate the contribution 

'^ii'^id = ^ La^ijLa^M {{ciCjCkCi, T^} - 2ciCjTjj:kCi) (A.3) 

ijkl 

where we used the hermiticity L^^ = of the Lindblad operators. The expression in 
the parenthesis can be rewritten as CiCj [ckCi, r,J — [cjCj, c^q and using the symmetry 
under exchanging the pair of indices (i, j) ^ {k, I) one finds 

^ii^u) = ^ X] La,ijLa,kl [CiCj, [CkCl, Tj] . (A.4) 
ijkl 

Finally, applying (1A.1I) twice and using L^^ki = —La,ik one arrives at 

= \Y^ L^^ijL^^kic\cjc\ciT^ . (A.5) 

ijkl 

The full time evolution operator (E]) is then obtained by normal ordering, summing over 
a and adding the coherent contribution in Eq. (lA.2p . 

Appendix B. Solution of the Bethe equation 

The quantization of the wavenumber 9 is determined by Eqs. ( IT8l) with indices I = 
and / = — 1. The system of these two equations for the amplitudes A and B in 
( |T9l) has a nontrivial solution only if the determinant of the coefficients vanishes. This 
condition leads to the following Bethe equation 

i/3 sin ^ [cos ^A^ - = sin ^A^ (B.l) 

where the parameters are defined as 

A sing/2 1 f jy .^Njiyi 



(3 = 2--::^, R{q) = 1 z^Vg-.v. ^^-iVg..v. _ 2) 
7 cosg 2 V / 

The second equation then fixes the ratio of the amplitudes B/A = 821- Note, that the 
form of R{q) is set solely by the boundary condition. 
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The solution of the Bethe equation requires some attention since it depends on the 
parity of N as well as qj. The problem is essentially the same as the one treated in 
[16] only the exact form of the parameter /3 differs. This, however, does not change the 
qualitative form of the spectrum and leads to the same families of eigenvalues. In the 
following we give a short summary of the results on the spectrum and refer to [16] for a 
detailed analysis. For the sake of concreteness we choose iV = 4/c + 2 with an arbitrary 
integer k. 



The complex eigenvalues Oj 

The first family is obtained by solving 



./!sm(-=.i 'ft (B,3) 

' tan ^ qj even 

and the scattering phase is given by 

^21 = -^, i = 13 Sine. (B.4) 

In general, Eq. (IB.3|) yields complex solutions for 9 which have different asymptotic 
expressions depending on the value of ^. For |^| < 1 one has = Oq^ + i(5< where the 
real and imaginary parts read 

_ / 5(2n+l), n = 0,l,...,f -1 . l,„l+«< 
n = l.....f-l ■ '<-n"'^—. 

with ^< = /3sin6'o<. The upper and the lower solutions refer again to qj being odd and 
even, respectively. For |^| > 1 the solutions 9^ = + are given by 

n / #2n, n = l,...,f-l 1 e> + l 

'°>^l|(2n + l), n^O,l,...!f-l ' '>-n''^I^I (^-^^ 

with = /3sin6'o>. Because of the condition |^| > 1, the solutions 6y only exist in the 
interval O^. < 9q < -k — O^. with O^. = arcsin In case 6c < tt/N, two special solutions 

for qj odd appear at 9 = 6 and 9 = n — 6 with 6"^ = -j^. Note, that the corrections to 
these asymptotic forms are 0{1/N'^) except from the vicinity of the branch cut |^| ~ 1 
where the approximation is not valid. 



The real eigenvalues 9jj 

Equation (]B.1|) can also be satisfied by setting cos N9 = R{q) = ±1 where the ± sign 
refers to qj being odd and even, respectively. The solutions are 

{—2n, n = 1, — — 1 

|(2n + l), n = 0;i,'.'..',f-l • (^-^^ 

The scattering phase is exactly 5*21 = —1 and therefore this family constitutes the 
antisymmetric solutions of the potential scattering problem. The latter property is also 
evident from the fact that the solutions are independent of /3 and thus of the exact form 
of the potential. 
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The diffusive eigenvalues 6^ 



There exists a special family of eigenvalues which in turn correspond to the bound state 
solutions of the Bethe equation. Setting 6*^ = ±| + i?7 one has 



cosh?7 



tanh ^^ Qj odd 
coth ^ qj even 



(B.8) 



where the sign ± refers to the cases /3 > and /3 < 0, respectively. In the N ^ oo limit 
the solutions are given by 

7] = acosh|/3|-i (B.9) 

^ 1, finite size corrections St] are exponentially small 



and, except from the vicinity | 
in N 

g— A''acosh 1/3 



(B.IO) 



Therefore ^ — )■ 1 and the scattering phase vanishes, corresponding to exponentially 
decaying, localized eigenmodes. 

The diffusive eigenvalues have to satisfy \f3\ < 1. Thus, for 7 < 2 A, they only exist 
in the intervals < g < gi and 27r — gi < q < 2tt where qi is determined by the condition 
j3{qi) = 1. For larger values of 7 one has an additional interval q2 < q < 2tt — q2 where 
P{Q2) = — 1- The above conditions are sketched in Fig. IB II for the same values of 7 as 
used in Fig. [21 see text. For 7 ^ A the wavenumber qi — )■ 7r/2 from below as well as 
q2 — )■ vr/2 from above and one has 6d eigenvalues almost everywhere in the spectrum, 
recovering purely diffusive behaviour. 



3 
2 
1 

P 




1^0.1 
7=3 





Figure Bl. Plots of the function /3 in Eq. (jB.2p for different values of 7 and A = 1. The 
solutions gi of /3 = 1 are indicated in corresponding color. For 7 = 3 the wavenumber 
(72 indicates the solution 13 — ~l. 
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